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Using 106 million (3686) events collected with the BESIII detector, we measure multipole amplitudes for 
the decay (3686) — Yxe1,2 + yyJ/w beyond the dominant electric-dipole amplitudes. The normalized 
magnetic-quadrupole amplitude for 1»(3686) > xc,» > yyJ/wv and the normalized electric-octupole am- 
plitudes for (3686) — YXec2, Xe2 — yJ/w are determined. The M2 amplitudes for (3686) > yxe1 
and Xc1,2 — yJ/w are found to differ significantly from zero and are consistent with theoretical pre- 
dictions. We also obtain the ratios of M2 contributions of (3686) and J/x) decays to Xe1,2, b3/b3 = 
1.35 + 0.72 and al/a2 = 0.617 + 0.083, which agree well with theoretical expectations. By consider- 
ing the multipole contributions of x«1,2, we measure the product branching fractions for the cascade decays 
(3686) + 7xc0,1,2 > yyJ/w and search for the process 7. (25) — yJ/w through (3686) — ^n«(25). 
The product branching fraction for (3686) — yx-0 — -yyJ/w is 3o larger than published measure- 
ments, while those of 7)(3686) — xc» — yyJ/w are consistent. No significant signal for the decay 
(3686) > yne(29) — yyJ/w is observed, and the upper limit of the product branching fraction at the 90% 


201711.02290v1 


chinaXiv 


confidence level is determined. 


PACS numbers: 14.40.Pq, 13.20.Gd, 13.40.Hq 


I. INTRODUCTION 


The processes (3686) — 71Xe1,2 and Xe1,2 — Y2I/W 
are dominated by electric-dipole (E1) amplitudes but allow 
for higher multipole amplitudes as well, such as the magnetic- 
quadrupole (M2) and electric-octupole (E3) transitions. The 
contributions of these higher multipole amplitudes give infor- 
mation on the anomalous magnetic moment « of the charm 
quark [ibl and on the admixture of S- and D-wave states (3h. 
The normalized M2 contributions for (3686) — y1 Xc1,2 and 
Xe1,2 — "yaJ/ , which are referred to as pi? and al? with 
the superscript representing x«1,», are predicted to be related 
to the mass of the charm quark, mec, and & n . By assum- 
ing m, — 1.5 GeV/c? and ignoring the mixing of S- and D- 
wave states, the contributions i^ and a. corrected to first 
order in Ey, ,/m., where E», , is the energy of 71,2 in the 
rest frame of the mother charmonium state, are predicted (4) 
to be 


b} = En [9 (3686) > vixe q + K) = 0.029(1 + &), 
Ame 

E c 
al = -Faka ct la) = —0.065(1 + k), E 
p= Em [&(3686) > Xixezl 4 | w) = 0.029(1 + x), 

V5 4me 

3 Ey [Xe2 > yaJ/b 

a2 = Z val T / la HK) = —0.096(1 + &), 


respectively. The ratio of the M2 contributions of 1/(3686) — 
1Xe1 to (3686) > 71Xc2 (xa — "aJ / to Xe2 > Y2T/V) 
is independent of the m, and « of the charm quark to first or- 
der in E,/m, and predicted to be 53/52. = 1.000 + 0.015 
and al/a2 = 0.676 + 0.071, respectively [5], where the 
dominant uncertainties come from ignoring contributions of 
higher-order in (E.,/m.)?. Higher order multipole ampli- 
tudes can be obtained by investigating the angular distribu- 
tions of the final-state particles iid . Several experiments 
have searched for higher-order multipole amplitudes ra 


[12]. The CLEO experiment reported significant M2 contribu- 


tions in (3686) > ^1xc1 and yc,» — Y2I/y by analyz- 
ing 24 million (3686) decays [5]. Recently, BESIII found 
evidence for the M2 contribution in (3686) — xc» with 
Xo ntr /K*K- (12). 

The experimentally observed charmonium states and their 
decay can be reproduced reasonably well by calculations 
based on a potential model and by perturbative quantum chro- 
modynamics [13]. However, for the El radiative transitions 
of 1»(3686) — 71Xc0,1,2, there are significant discrepancies 
between different model predictions Hlg and the Parti- 
cle Data Group (PDG) average (17). The partial widths of 
w(3686) + 71Xc0,1,2 are predicted to be 26, 29, and 24 keV, 
respectively, by using the Godfrey-Isgur model {id}, which 
deviate by — (13 + 3.5)%, (1.4+ 4.6)96, and —(11.8+ 3.9)% 
from the averages of experimental measurements (171. 

In this paper, we report on a measurement of the higher- 
order multipole amplitudes in the processes of 1 (3686) — 


Y1Xc1,2;Xc1,2 > Y2I/W, where the J/z is reconstructed in 
its decay modes J/y — 4+4 (€ = e/p). The measurements 
make use of the joint distributions of the five helicity angles in 
the final-state. Using the invariant mass of y2J/w, we obtain 
the product branching fractions of (3686) — ¥1Xc0,1,2 > 
^1^yaJ /' and search for (25) > *3J/1» produced through 
(3686) — «19.(25). In the measurement of the product 
branching fractions of 7(3686) > Y1Xc0,1,2 — o*9J/w. 
the multipole contributions of x,» are considered for the first 
time. The results presented in this manuscript supersede the 
ones in Ref. [18]. The analyses are based on a sample of 
156 pb^! taken at a center-of-mass energy 3.686 GeV, cor- 
responding to 106 million (3686) [19]. A 928 pb^! data 
sample taken at 3.773 GeV and a 44 pb~! data sample 
taken at 3.65 GeV are used to estimate the backgrounds from 
QCD processes. 


Il. BESIII DETECTOR AND MONTE CARLO 
SIMULATION 


The BESIII detector is described in detail in Ref. Rih. It 
is an approximately cylindrically symmetric detector which 
covers 93% of the solid angle around the collision point. The 
detector consists of four main components: (a) a 43-layer 
main drift chamber provides a momentum resolution of 0.5% 
for charged tracks at 1 GeV/c in a 1 T magnetic field; (b) 
a time-of-flight system (TOF) is constructed of plastic scin- 
tillators with a time resolution of 80 ps (110 ps) in the bar- 
rel (end caps); (c) a 6240 cell CsI(TI) crystal electromagnetic 
calorimeter (EMC) provides an energy resolution for photons 
of 3.096 (5.096) around 0.3 GeV in the barrel (end caps) [23]: 
(d) a muon counter consisting of nine (eight) layers of resis- 
tive plate chambers in the barrel (end caps) within the return 
yoke of the magnet with a position resolution of 2 cm pro- 
vides muon/pion separation. A GEANT4 based detector 
simulation package has been developed to model the detector 
response used in Monte Carlo (MC) generated events. 

A MC simulated sample of 106 million generic (3686) 
decays ("inclusive MC") is used for general background stud- 
ies. The (3686) resonances are produced by the event 
generator orc fi The known decays are generated by 
BESEVTGEN [25] with branching fractions taken from the 
PDG Ii while the remaining decays are generated accord- 
ing to the LUNDCHARM model . Exclusive MC sam- 
ples for signal decays are generated to optimize the selec- 
tion criteria and to determine the detection efficiencies. The 
(3686) > 7xXc0,1,2 — YyyJ/1 decays are generated with 
angular distributions determined from data, and the ne(25) > 
yJ/w decay is generated according to the HELAMP model 
in EVTGEN i25]]. To estimate the background contributions 
from 1/(3686) decays, the exclusive MC samples (3686) > 
1J [Ab, m9 J [ap, n?n? J [Ab yvy Ji are generated according to 
the HELAMP, JPIPI 251]. and PHSP models, respectively. 
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To investigate QED processes backgrounds, radiative Bhabha 
and dimuon events (e*e^ — e*e /p* p) simulated with 
BABAYAGA V3.5 A as well as (3770) > ¥ and 
^"rsn (3686) — yx«;, 19 J/V» produced by KKMC , are 
used together with the experimental data at 3.773 GeV. 


IIl. EVENT SELECTION 


The signal decay (3686) —  ^ixeoi(m«(29)) — 
^n*yaJ[ib, Jf — Cre (€ = e, p) consists of two charged 
tracks and two photons. Events with exactly two oppositely 
charged tracks and from two up to four photon candidates are 
selected. Charged tracks are required to originate from the 
run-dependent interaction point within 1 cm in the direction 
perpendicular to and within +10 cm along the beam axis and 
should lie within the polar angular region of | cos0| < 0.93. 
The momentum p of each track must be larger than 1 GeV/c. 
The energy deposit E in the EMC and E/p of each track are 
used to identify muon or electron candidates. Tracks with 
E < 0.4 GeV are taken as muons, and those with E/p > 0.8c 
are identified as electrons. Events with both tracks identi- 
fied as muons or electrons are accepted for further analy- 
sis. Photons are reconstructed from isolated showers in the 
EMC, where the angle between the positions in the EMC of 
the photon and the closest charged track is required to be 
larger than 10 deg. The energy deposited in the EMC is cor- 
rected by the energy loss in nearby TOF counters to improve 
the reconstruction efficiency and the energy resolution. The 
energy of each photon shower is required to be larger than 
25 MeV. The shower timing information is required to be in 
coincidence with the event start time with a requirement of 
0 < t < 700 ns to suppress electronic noise and showers 
unrelated to the event. 

A four-constraint (4C) kinematic fit is performed for the 
two lepton candidates and all possible two photon combina- 
tions with the initial 1»(3686) 4-momentum as a constraint. 
If more than one combination is found in one event, the one 
with the smallest x2c value is kept. The x2c is required to 
be x2c < 60, where the requirement is determined by opti- 
mizing the statistical significance S/V S + B for the n. (25) 
channel. Here, S is the number of events in the 7, (2) signal 
region 3.60 < M^C€(454* £7) < 3.66 GeV/c? (y2 denotes the 
photon with larger energy, and M^€ is the invariant mass with 
the energies and momenta updated with the 4C kinematic fit) 
obtained from the exclusive MC sample, and B is the num- 
ber of corresponding background events determined from the 
106 million inclusive MC sample and a continuum data sam- 
ple collected at a center-of-mass energy of 3.65 GeV. The lat- 
ter is normalized to the luminosity of the (3686) data sam- 
ple. The branching fraction of the decay 7.(2S) > y2J/w is 
assumed to be 1%. 

To select events including the J/ọ% intermediate state, the 
invariant mass of the lepton pair is required to be in the re- 
gion of 3.08 < M*€(/*£7) < 3.12 GeV/c?. In addition, 
to remove (3686) > 7? J/y and (3686) — J/w back- 
grounds, events with an invariant mass of the photon pair in 
the regions 0.11 < MC (yy) < 0.15 GeV/c? or M*€ (qw) > 


—- data 
B Background 
—— X, Signal MC 


oe’ 
t 
eee 


Events / 2 MeV/c? 
a, 


3.46 3.48 3.5 3.52 3.54 3.56 3.58 3.6 
M'^(y I) (MeV/c?) 


FIG. 1. Mass distributions of M^*€ (yalt) for events in the x«1,2 
region. Black dots correspond to data, and red histograms are ob- 
tained from the signal MC samples scaled by the maximum bin. 
The green dashed histogram is the background contribution obtained 
from the inclusive MC samples. The arrows denote the signal re- 
gions. 


0.51 GeV/c? are rejected. A MC study shows that this re- 
moves 97.9% of the 1? J/1» events and almost 100% of the 
nJ/w events, while the efficiencies of the signal channels 
for xco, Xci, Xe2, and (25) are 74.7%, 90.096, 93.9%, and 
88.0%, respectively. 


IV. MEASUREMENT OF HIGHER-ORDER MULTIPOLE 


AMPLITUDES 
Figure shows the M#°(y2¢+é-) invariant- 
mass distribution for the selected x1, candidates. 
The signal regions for x. and Xx.» are defined 


as 3.4.96 < — M'C(wSt*07) — «3.533 GeV/c? and 
3.543 < M*C(yoltl-) < 3.575 GeV/c?, respectively. 
We find 163922 x., candidates and 89409 x.5 candidates. 
The background is estimated from the inclusive MC sample. 
The total number of background events is found to be 1016 
(0.7%) within the x, signal region and 883 (1.0%) in the x.» 
region. For the x«i (xc2) channel, the dominant background 
is the contamination from xco (xc1). Some backgrounds 
stem from (3686) — yyJ/w and r°r°J/w decays. The 
QED process ete~ — (*-^jsn ran contributes about 109 
events for Xcı and 135 events for Xe2. Non-J/ background 
is negligibly small according to the sideband analysis. 


Events in the signal regions are used to determine 
the higher-order multipole amplitudes in the (3686) — 
^iXei2 — ^yaJ/* radiative transitions. The normalized 
M2 contributions for the channels (3686) — ^1x«1,» and 
Xe1,2 > y2I/ are denoted as m and as", respectively. In 
the x.» decays, the E3 transition is also allowed. The corre- 
sponding normalized E3 amplitudes are indicated as 52 and a2 
for (3686) > ^1xc» and xeo — *2J/ v, respectively. 
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IVA. Fit method 


We perform an unbinned maximum likelihood fit to obtain 
the higher-order multipole amplitudes following the proce- 
dure as described in Ref. [12]. The log-likelihood function 
is built as In £L; = ln £ — In £j, where £L = Mi E sali) 
denotes the product of probability densities for all candidates 
in the signal region, N is the number of the candidates, and F 
is the probability density functions (PDFs). The contribution 
to the likelihood from background events, Lp, is estimated us- 
ing the inclusive MC sample and continuum data. 

The PDFs F for the joint angular distribu- 
tions of the x1» decay sequences are defined as 
Wy... (01,02,02,05,03,a3 3,03 3) 

Ws, 7 (42,312 3) f 
Wy. (01, 02, 2,03, bs, 02 3; bj), is derived from the he- 
licity amplitudes and the Clebsch-Gordan relation ih, while 
Wys (a3 3,b3 3) is used for the normalization. 6; is the 
polar angle of ^ in the (3686) rest frame with the z axis 
in the electron-beam direction. 05 and $» are the polar and 
azimuthal angles of 72 in the x; rest frame with the z axis 
in the ^4 direction and ¢2 = 0 in the electron-beam direction. 
05 and #3 are the polar and azimuthal angles of £* from 
J/w —5 C C7 in the J/v rest frame with the z axis aligned to 
the y2 direction and $5 = 0 in the 7 direction. 

The formula Wy., (01, 02, 62, 03, $3, az 3, bd 4) for the he- 
licity amplitudes has been discussed in Refs. [28]. 
Using the same method as reported in Refs. 15}, [u 12, the 
joint angular distributions W,.., can be expressed in terms of 
až 3 and b7 4 as 


The term in the numerator, 


Wx.; (04, 05, Q2, 03, $3, 03 3; b; 3) 


=D ni AR Bion Bi = 
with 

Aj\ _ (V0.5 V0.5 at 

Al] \V05 -V05/) (a3) 

BjV _ (V0.5 v0.5 (bi 

Bi) ~ \ vos —V05) NJ 

A2 /3 ^15 2V3 a? 

9 1 1 (3) 

A? | =—=| 3 vb 4 az | , 

Aj] V9 Asy2 —V10 v2) Nai 

B2 1 [v3 v5 2V9V fv 

B? | = — | 3 v5 —4 i 

Bj V9 Asyg -v0 v2) S 
where By, and Bi are the helicity amplitudes for 


(3686) — mx. Aj, and A7 are those for Xey > 


|v| 
y2I/p. y (a1? + (a5)? = 1, y (a1)? + (a2)? + (a3)? = 1, 
and similarly for bis. The coefficients an (27-3 fer xa. ) 


—1,...,96 for xc» 
are functions of 01,05, 2,03, 03. For the normalization, 


high-statistics phase-space (PHSP) MC samples are gener- 
ated. 
The normalization factor is expressed as 


Wy.; (a3 3, by 4) 


XR Wes (01 (i), 02), doi), 05i), $2 (1), a2 5, 63.3) 
Np 


=> TrA Ap Bp Bin 
(4) 


where N p is the number of selected events. In such a way, the 
detector efficiency is considered in the normalization. 


IV.B. Fit results 


By minimizing — In £+, the best estimates of the high-order 
multipole amplitudes can be obtained. To validate the fit pro- 
cedure, checks are performed with MC samples for x1,» sep- 
arately, where the MC samples are generated based on a pure 


El transition model (nois — 0, Me = 0) or an arbitrary 


higher-order multipole amplitude (25/3 z 0, n z 0). The 
fit values are consistent with the input values within 1c of sta- 
tistical uncertainty. An unbinned maximum likelihood fit to 
the joint angular distribution for data is performed, and the 
corresponding angular distributions are depicted in Fig. []to- 
gether with the relative residual spectra. The fit results are 
listed in Table [I| where the first uncertainties are statistical 
and the second ones are systematical as described in Sec. [VI] 

The statistical significance of a nonpure E1 transition is cal- 
culated to be 24.50 (13.50) for xc (x«2) by taking the dif- 
ference of the log-likelihood values for the fits with higher- 
order multipole amplitudes included and fits based on a pure 
EI transition, taking the change in the number of degrees 
of freedom, Andf = 2 (4), into consideration. Similarly, 
the statistical significance of the E3 contribution for Xe2 is 
2.30, as obtained by comparing the log-likelihood values 
between the nominal fit and a fit based on the assumption 
that E3 contribution is zero. A Pearson-X? test is per- 
formed to validate the fit result. Each angular dimension 
(1.e., cos 01, cos 05, b2, cos 05, $3) is divided equally into eight 
bins. This leads to a total of 8° = 32768 cells. The X? is de- 
fined as 


ayn oP Ee ca 
x= nDT + pBKG , 
where ppt is the number of events in the ith cell for data, 
nBKG is the number of the background contribution deter- 
mined by the inclusive MC sample, and nM€ is the number 
of events for the luminosity-normalized MC sample produced 
according to the best fit values for aj 3 and bj 3. The number 
of events of the MC sample is 40 times larger than of the data. 
For cells with fewer than ten events, events in adjacent bins are 
combined. The test results in x? /ndf = 9714.7/9563 = 1.02 
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TABLE I. Fit results for als and 57 5 for the process of (3686) > 71Xc1,2 > Y1Y%2J/%; the first uncertainty is statistical, and the second is 


systematic. The pi. 303.5 


for x4 and x? /ndf = 5985.2/5840 = 1.02 for x3, demon- 
strating that the fit gives an excellent representation of the 
data. 


V. MEASUREMENT OF 
B(V(3686) — yxcs — vyJ/V) AND SEARCH FOR THE 
PROCESS n«(28) — yI /4 


With the selected e*e^ — ^72J/1^ candidates, we mea- 
sure the product branching fractions of the decay (3686) — 
^nXc0,1,2 — ^nyaJ/ and search for the process rg. (29) > 
y2J/p. For the J/y — ete~ channel, additional require- 
ments are applied to suppress the background from radiative 
Bhabha events [e*e^ = YISR/FSRE ET, where YISR/FSR 
denotes the initial-/final-state radiative (ISR/FSR) photon(s)]. 
Since the electron (positron) from radiative Bhabha tends to 
have a polar angle cos0,.(,-) close to +1 (-1), we apply a 
requirement of cos @,+ < 0.3 and cos ĝe- > —0.3. These re- 
quirements suppress 77% of the Bhabha events with a reduc- 
tion of the signal efficiency by one-third. The corresponding 
MC-determined efficiencies are listed in Table[IT] 

A 4C kinematic fit has the defect that the energy of a fake 
and soft photon will be modified according to the topology 
of a signal event due to relatively large uncertainty, which re- 
sults in a peaking background signature in the M^€ (^5 J/») 
invariant-mass spectrum. To remove the peaking background, 
such as radiative Bhabha and radiative dimuon (e*e^ — 
YISR /FSRAT pi), a three-constraint (3C) kinematic fit is ap- 
plied, in which the energy of the soft photon (71) is left free in 
the fit. The detailed MC studies indicate that the 3C kinematic 
fit does not change the peak position of the invariant mass for 
signals and the corresponding resolutions are similar to those 
with the 4C kinematic fit. 


V.A. Background study 


The backgrounds mainly come from v/(3686) transitions to 
Jib and from e*e^ — (UU nsnrsn (€ = e/ nu). The 
other background, including (3686) > nJ/w, ^en J / 1» and 
non-J/1» backgrounds, is only 0.3% of that from w(3686), 
which is neglected. 

The backgrounds from (3686) transitions to J/«» in- 
clude (3686) + yy J/p, n?n? S/W, n9 J/v. High-statistics 
MC samples of these decays are generated to determine 


are the correlation coefficients between a33 and b33. 
al = —0.0740 + 0.0033 £ 0.0034, 63 = 0.0229 + 0.0039 £ 0.0027 
Phab, = 0.133 
az = —0.120 £ 0.013 £ 0.004, 63 = 0.017 Æ 0.008 £ 0.002 
az = —0.013 + 0.009 + 0.004, b2 = —0.014 + 0.007 + 0.004 
Pagta = —0.605, Pasas = 0.733, Pandy 
Pagby = —0.422, poop, = 0.384, P2303 = —0.024 


= —0.095 


their distributions and contributions. With the published 
branching fractions (17h, which have been measured precisely 
by different experiments, the estimated number of events 
for (3686) — 1°n°J/w,2°J/w and the efficiency for 
w(3686) — yyJ/w are obtained as summarized in Table [I] 

The second major source of background  in- 
cludes radiative Bhabha and dimuon processes, 


ete” — LHL ysR/Frsr(VIsR/FsR) and (3686) — 
+l- tysn(orsn) (l= e/u). To precisely describe the shape, 
the background is divided up into two parts: (* ^ with one 
radiative photon and ¢+£~ with two radiative photons. For 
the background from 7)(3686) > £* (^ ^rsn (^n), the ratio 
of event yields between the two parts (Ne+¢—~y/Ne+e-) is 
obtained by a MC simulation. For the background from ra- 
diative Bhabha/dimuon processes, the ratio Ny+¢—+-/Ne+e- 
is obtained by a fit to a 928 pb~! data sample taken at a 
center-of-mass energy of 3.773 GeV. After the event selection 
imposed on the data, the remaining events are mainly radiative 
Bhabha/dimuon events, and a small contribution originates 
from (3770) + yxez and decays of (3686) produced in 
the ISR process. In the fit, the shapes of the M°° (72+ £-) 
distributions for the Bhabha/dimuon processes are determined 
from a 1/(3686) > £* U^ ygsn (yrsn) MC sample by shifting 
the M3C (42427) from 1/(3686) to (3770) according to 
the formula m’ = a * (m — mo) + mo, where mo = 3.097 
GeV/c? is the mass threshold of y.J/w, and the coefficient 
a = (3.773—mpo)/(3.686—mo) = 1.15 shifts the events from 
3.686 to 3.773 GeV. The shapes of the backgrounds are based 
on MC simulation, while the amplitude of each component 
is set as a free parameter. Thus, the cross section weighted 
ratio of the backgrounds ete~ — €*é~~ygR/rsr(V1SR/FSR) 
and (3686) — LL- psr(yrsr) for the two parts is 
Nete=ya/Nete=y = 1.2038 00081 (IN, e ui s IN uu. = 
0.689 + 0.044) for the e*e^ (ut) channel. The quan- 
titative results and shapes will be used in the simultaneous 
fit. 


cS 


V.B. Simultaneous fit to M?S (y2e £^) 


Figure [3] shows the M°°(72¢+¢-) distributions for se- 
lected candidates of the two channels of J/y — ete~ and 
J/v — wtp, where clear signals of Xco,1,2 can be ob- 
served. No evident 7.(25) signature is found. A simultane- 
ous unbinned maximum likelihood fit is performed to obtain 


Data 


—Fit MC 
- pure E1 


- Data 


—Fit MC 
- pure E1 


--Data 
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- pure E1 
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FIG. 2. Results of the multidimensional fit on the joint angular distribution and the projections on cos 01, cos 0», cos 05, $2, 63 of the final- 
state particles. The upper ten plots show the angular distributions for the x«1 channel, and the lower ones are for the y-2 channel. The black 
dots with error bars represent data subtracted by background, the red histograms are the fit results, and the blue dashed lines are pure E1 
distributions. The lower plots depict the relative residual x = (Naata — Ng«)/ V/ Naata of the fit. 


TABLE II. Detection efficiencies (c) for channels of (3686) > ^/xo,1,2, 1c (28), yyJ/w and the number (N) of estimated background for 


channels (3686) — m? J/y, n^n? J/w scaled by the decay branching fraction and the total (3686) number. 
Channel exco (96) exer (96) Exea (96) eg. (23(V6) Errare (%) IN o; 


, N 


ab nOn T ab 


ee” 15.1 
uu 32.7 


20.1 
44.1 


20.3 
44.0 


the signal yields. The common parameter for the two J/w de- 
cay channels is the product branching fraction (BProduct) of 
the cascade decays (3686) > 7xXc0,1,2(Ne(2S)) > yyJ/w. 
The number of signal events for each channel is NY(°89) x 
Brroduct x B(J/y —> (* 0-7) x e. In the fit, the branching 
fractions for J/v — e*e™/u*u and the total number of 
~(3686) events are fixed to the values in Refs. and [19], 
respectively. The efficiency e is obtained from the signal MC 
sample with the higher-order multipole amplitudes considered 
as listed in Table[I]] The fit contains three Xc0,1,2 components, 
the 1. (25), and the background. The signal line shapes of the 
Xc0,1,2 are parametrized as 


(6) 


where BW (m) is the Breit-Wigner function for Xco,1,2 with 
the masses and widths fixed at their world average values ei 
R represents the mass resolution, and e(m) is the mass- 
dependent efficiency. The product [BW (m) & R x e(m)] 


(ES, x E3, x (BW(m) & R x e(m))) ® G(u, o), 


16.9 
37.0 


17.1 
38.0 


26.8+0.7 246.5+4.5 
65.2+1.7 500.9+9.1 


can be directly determined from the MC simulation, where 
the MC events are generated with the simple Breit-Wigner 
function using the higher-order multipole amplitudes with the 
angular distributions of the final-state particles. £54 is the 
energy of the radiative photon 71 of (3686) > 1x. in 
the (3686) rest frame, and E» is the energy of the y2 of 
XcJ > aJ [b in the x; rest frame. The factor E315 stems 
from the two-body PHSP and the EI-transition factor, and the 
Breit- Wigner function modified by the E215 factor is for the 
XcJ invariant-mass distribution. The line shape is convoluted 
with a Gaussian function (denoted as G) accounting for dif- 
ferences in the invariant mass and mass resolution between 
the data and the MC simulation. The mean u and standard 
deviation c of the Gaussian functions are obtained from the 
fit to the data in a region of [3.36 < M®°(y20tl-) < 3.61 
GeV/c?] by assuming no dependence between the e*e^ and 
pt decay modes as well as between Xco,1,2. The results 
indicate jj < 0.35 MeV/c? and ø € 0.73 MeV/c?. Similarly, 
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the signal line shape of the ne (25) is described by 


(ES, x BY, x (B(m) & R x €(m))) @G(u,0), (D) 
where ER, represents the two-body PHSP and the MI- 
transition factor for (3686) — ^1n.(2S) and E; is the 
two-body PHSP and hindered M1 transition factor ‘itd Bo Bo 
for 7-(2S) + yaJ/ v. The [B(m) & R x e(m)] is also deter- 
mined by MC simulation with Es mass and width of n.(25) 
set to the world average values [17]. Since the mass of 7. (25) 
is close to those of Xey, the u irm o of the Gaussian are fixed 
to the values obtained from a fit to the Xco,1,2 signals only. 

The shapes of backgrounds | (3686) — 
19 J /b, n9 n9 Jb, yyJf and e*e-(—2 — w(3686) — 
LELT msnrsn (visn/FsR) are taken from MC sim- 
ulations. The numbers of ~(3686) — 7°J/w and 
(3686) —>  nx9«9J/w qr are fixed to the expec- 
tations as given in Table [I| For the background from 
e*e (2 w(3686) > M. the ratios of 
Ne+e-yy/Ne+e-y are fixed to 1.203 for the e*e^ chan- 
nel and to 0.689 for the ~*~ channel as described 
above. In the fit for the final results in the region 
(3.36 < M?C€(q4£*07) < 3.71 GeV/c?), the parame- 
ters of the smearing Gaussians for Xc0,1,2 and (25) are 
fixed, while the numbers of events for ¥¢0,1,2 and (25), 
v(3686) > yyJ/y, ete” — *U msnyrsn(msnyrsn) 
are free parameters. Figure [3] shows the M?C(yo£*£-) 
distributions, the results of the unbinned maximum likelihood 
fit, and the relative residuals. The x? /ndf of the fit is 1.88 for 
the u+ u~ channel and 1.83 for the e* e^ channel. 

The product branching fractions from the fit are (15.8 + 
0.3) x 1074, (351.8 + 1.0) x 1074 and (199.6 + 0.8) x 1074 
for Xco,1,2 With statistical uncertainty only, respectively. The 
branching fraction of (3686) — yyJ/w is determined to be 
(3.20.6) x 1074. All measured branching fractions are con- 
sistent with the previous measurement of BESIII [18 [18]. Since 
no significant 7.(2S') signal is found, an upper limit at the 
90% C.L. on the product branching fraction is determined by 
a Bayesian approach using a uniform prior, i.e., finding the 
values corresponding to 90% of the probability distribution in 
the positive domain. 


VI. SYSTEMATIC UNCERTAINTIES 


The main sources of systematic uncertainty for the mea- 
surements of higher-order multipole amplitudes are the un- 
certainties in the efficiency, the kinematic fit procedure, the 
fit procedure of the combined angular distributions, statistical 
fluctuations of the MC sample, and the background contami- 
nation. 

A simulated sample of events distributed uniformly in 
PHSP is used to normalize the function Wy, ,. A differ- 
ence of detection efficiencies between the MC sample and 
the data will result in a shift in the measurement, which is 
taken as the systematic uncertainty. From the studies of the 
tracking efficiency for electrons and muons with the control 
samples of (3686) 3 ntr J/wv,J/v — ete /n*u- 


decays, and the photon efficiency with the control samples 
from 7)(3686) — 2(r*«-)n? decays and radiative dimuon 
events, the difference in the detection efficiencies between the 
data and MC is found to be polar angle dependent with the 
largest value 0.006 + 0.003, which may change the helicity 
angular distribution. The corresponding effect on the higher- 
order multipole measurement is estimated by varying the effi- 
ciency with an asymmetric function of cos 0,.. and cos 0.,, as 
p(cos 0.1, cos 05.) = (1.0+0.003 cos 6,1 —0.006 cos? 0,1) x 
(1.0 + 0.003 cos 62+ — 0.006 cos? 6,4) [which corresponds to 
a 0.9% (0.3%) difference for cos 0 = —1 (1); 6,4 is the polar 
angle for one photon, and 0,.- is for one charged track]. Twice 
the difference with respect to the nominal result is taken as a 
systematic uncertainty. For the kinematic fit, the track helix 
parameters are corrected to reduce the difference in the x2c 
distribution between the data and the Eu simulation accord- 
ing to the procedure described in Refs. (31,82). These PHSP 
MC samples without and with the helix a. are used to 
normalize Wx. ,, respectively, and the resultant difference is 
taken as the systematic uncertainty. 

To estimate the uncertainty from the fit procedure, 200 MC 
samples using the high-order multipole amplitudes are gen- 
erated, followed by a complete detector simulation. Each 
sample has 165 thousand (90 thousand) selected events 
for Xcı(Xc2), and the same multipole analysis procedure 
is applied for each sample. The differences in al, bl 
(a3, a3, b2, b3) between the input and fitted values are 
Gaussian distributed. The mean values of the Gaussians 


are Ji = (2 + 3) x 1074, Hos = (—6 + 3) x 10-4 
(Haz = (17+ 13) x104, Haz = (—4 8) x 1074, Hy = 
(16 + 6) x 1075, m3 = (—32 + 7) x 1074) and are taken 


as the systematic uncertainty. The statistics of the MC sample 

for the normalization, about 3.6 (1.8) million events, may af- 

fect the fit results. For the normalization function, Eq.(4). the 
B — Han n 9 for xci : 

variance for an (5,36 for 3.4) IS 


Vin) = (C - eget. 

The standard deviation for each coefficient is o(@n) = 
V V(a;). The largest change in parameters aj and bl by 
varying the coefficient by +1o for the x«i channel (a2, a3, 
and b3, 02 for the x.» channel) is taken as the systematic un- 
certainty. 

The main backgrounds for the x«i channel come from 
(3686) — "xc, YXe2, n9 n9 J b, yyJ/Y, which contribute 
about 0.7% of the candidates according to a MC study. For the 
Xc2 channel, the main backgrounds come from 10(3686) — 
"YXcos VXe1, TOTP J [Ab yyJ/%, and the contribution is about 
1%. In the nominal fit, the contribution of background is esti- 
mated by the inclusive MC samples. To estimate the system- 
atic uncertainty, high-statistics MC samples for backgrounds 
are generated to redetermine the sha E and the contribution 
according to previous measurements [17 [18 Ba. The differ- 
ence in the fit results is taken as the systematic uncertainty. All 
the systematic uncertainties are summarized in Table [MI] The 
total systematic uncertainties are calculated by adding the in- 
dividual values in quadrature, thereby assuming that they are 
independent. 
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O FIG. 3. The results of a simultaneous maximum likelihood fit (top) and corresponding relative residual (Naata — Nnt )/V Naata (bottom). The 
© left panel is for the e* e^ channel, while right one is for the ut u~ channel. The black dots are the data, the blue curves are the fit results, 
A J and the red long-dashed lines are for y-0, 1,2 signals. The gray dashed, orange dot-dashed, and pink dotted lines are for backgrounds of 
A j (3686) — ^J /w, T? J/v, and n?n? J/1), respectively. The light-blue dot-dot-dashed and green dot-long-dashed lines are for backgrounds 


c with final-state particles composed of £* (^ y and £^ (^ yy. 


= TABLE III. The different sources of systematic uncertainties for the measurement of higher-order multipole amplitudes for the x1,» channels. 
Source EOE "y ) 
Efficiency of PHSP MC 14 
Kinematic fit 
Fitting procedure 
Statistics of PHSP MC 
Background 
Total 
IC The systematic uncertainties of the branching fractions and MC sample is 0.6% for the e* e^ channel and 0.1% for 
o measurement include uncertainties from the number of the yw" channel. To be conservative, we take 0.6% as the 


15 (3686) events (0.9%) Ed, the tracking efficiency (0.1% per 
lepton) [34], the photon detection efficiency (1.0% per pho- 
ton) [35], the kinematic fit, the J/1» mass window, the other 
selection criteria (N., < 4, veto 7° and rj, particle identifi- 
cation, cos@.+ < 0.3&&cos0,- > —0.3), the branching 
fraction of J/v — e*e- /u*q- (0.6%) i. the interfer- 
ence between 1(3686) + x«o — yyJ/w and nonresonant 
10(3686) — yyJ/wv processes, and the fitting procedure. 


The uncertainty from the kinematic fit is estimated by the 
same procedure as described in the multipole amplitude mea- 
surements. To estimate the uncertainty caused by the J/w 
mass requirement, a control sample in the x.1,» region 3.49 < 
M*C€(40*£7) < 3.58 GeV/c? is used. For data, the only 
background is from (3686) — 959 J/w, which is deter- 
mined in fitting with the exclusive MC shape. The efficiency 
of selection M^C(£*47) €(3.08,3.12) GeV is evaluated by 
comparing the number of signal events before and after the re- 
quirement, and the corresponding difference between the data 


systematic uncertainty. With the same sample, the systematic 
uncertainties related to the selection criteria Ny < 4, T°’ veto, 
7 veto, and leptons identification are also determined. The 
overall difference in the efficiency between the data and MC 
sample for these criteria is 1.6% and is taken as a systematic 
uncertainty. The additional systematic uncertainty due to the 
polar angle selection for the e*e^ channel is determined by 
varying the selection with +0.05 and fitting simultaneously 
again. The largest changes on the fit results are taken as the 
systematic uncertainty. 


To estimate the possible uncertainty from the interference 
between 1/(3686) — yyJ/ and v(3686) > yx. — 
yyJ/w, we repeat the simultaneous fit, taking the interfer- 
ence into account. The interference phase is found to be 
1.58 + 0.05. The changes in the signal yields are taken as 
the systematic uncertainty. Since the signal shapes are de- 
termined from MC simulation, the corresponding systematic 
uncertainty is estimated by an alternative fit with varying the 
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mass and width of y.0,1,2 with +1o of the world average val- 
ues for the signal MC shape. To estimate the uncertainty 
due to the background of (3686) > n? J/w, n?n? J/a and 
the ratio of Ny+¢-4-,/Ne+e-y for Bhabha and dimuon back- 
grounds, alternative fits are performed in which the numbers 
of expected background events (see Table [II] and the ratio of 
NA peg- / Nogea- are varied by +1ø. For Xco,1,2, the largest 
differences in the signal yields from the nominal values are 
taken as the systematic uncertainty. For the 7.(2.5') case, to 
be conservative, the one corresponding to the largest upper 
limit is taken as the final result. All systematic uncertainties 
of the different sources are summarized in Table[IV] The total 
systematic uncertainties are obtained by adding the individ- 
ual ones in quadrature, thereby assuming all these sources are 
independent. 


TABLE IV. Summary of all systematic uncertainties for the branch- 
ing fractions measurement. 


Source Xco (96) Xe1 (96) xca (95) me (2) (96) 
Nw (3686) 0.9 0.9 0.9 0.9 
Tracking efficiency 0.2 0.2 0.2 0.2 
Photon detection 2.0 2.0 2.0 2.0 
Kinematic fit 0.6 0.5 0.5 0.4 
J/w mass window 0.6 0.6 0.6 0.6 
Other selection 2.4 2:2 2.3 2.4 
B(Jlv > ete /u*u ) 0.6 0.6 0.6 0.6 
Interference 0.7 - - - 
Signal shape 0.7 0.9 1.0 - 
Background 0.1 0.1 0.1 - 
Total 3.6 3.4 3.5 3.4 


VII. RESULT AND SUMMARY 


Based on 106 million ~(3686) decays, we measure the 
higher-order multipole amplitudes for the decays (3686) > 
V1Xe1,2 — W172, /w channels. The statistical significance of 
nonpure E1 transition is 24.30 and 13.40 for the x. and xc 
channels, respectively. The normalized M2 contribution for 
Xc1,2 and the normalized E3 contributions for x.» are listed in 
Table[I] Figure[4] shows a comparison of our results with pre- 
viously published measurements and with theoretical predic- 
tions with m, = 1.5 GeV/c? and « = 0. The results are con- 
sistent with and more precise than those obtained by CLEO- 
c and confirm theoretical predictions [1 by. The M2 con- 
tributions for (3686) > ^ix« (bl), xai > Y2I/W (al), 
and xc3 — y2J/w (a3) are found to be significantly nonzero. 
The ratios of M2 contributions of Xc1 to Xe2 are independent 
of the mass me and the anomalous magnetic moment & of the 
charm quark at leading order in E-, /me. They are determined 
to be 


b3 /b2 = 1.35 + 0.72, 
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a} /a2 = 0.617 + 0.083. A 

The corresponding theory predictions are (b4/b3)in = 
1.000 + 0.015 and (a4/a%)tn = 0.676 + 0.071 (sl. By using 
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the most precise measurement of the M2 amplitudes a and 
by taking m, = 1.5 + 0.3 GeV/c?, the anomalous magnetic 
moment « can be obtained from Eq. (1). 


i4: Ame 1 

K = — Se 
E xa > y»J/V] ? (9) 
—1.140 + 0.051 + 0.053 + 0.229, 


where the first uncertainty is statistical, the second uncertainty 
is systematic, and the third uncertainty is from m, = 1.5 + 
0.3 GeV/c?. 

Based on the multipole analysis, we measure the product 
branching fractions for 1/(3686) > YXc0,1,2 > yyJ/w to be 
(15.8 + 0.3 0.6) x 1074, (351.8 + 1.0 + 12.0) x 1074, 
and (199.6 + 0.8 + 7.0) x 1074, respectively, where the 
first uncertainty is statistical and the second is systematic. In 
Fig. [5| the product branching fractions are compared to pre- 
vious results from BESIII , CLEO Bal, and the world 
average (17). The world average refers to the product of the 
average branching fraction of w(3686) — ^1x.; and the av- 
erage branching fraction of x«; — 72J/w, where the results 
of BESIII and CLEO are not included in the world average 
values. For all Xej, our results exceed the precision of the 
previous measurements. Compared to the previous BESIII re- 
sult, the results are consistent within 1c, but we have consid- 
ered the higher-order multipole amplitudes and improved the 
systematic uncertainty due to a more precise measurement of 
the total number of produced 7(3686) [19]. In addition, our 
measurement for the Xco channel is 3c larger than the result 
from CLEO and 3c larger than the world average value, while 
for the Xc1,2, our results are consistent with previous mea- 
surements. There are theoretical predictions for the branch- 
ing fraction (3686) — YXco,1,2 by several different mod- 
els without consideration of higher-order multipole 
amplitudes, which agree with each other poorly. The results 
in this measurement will provide a guidance for the theoretical 
calculations. 

We also search for the decay 7.(2S) — yJ/w through 
(3686) — yn-(2S). No statistically significant signal is 
observed. Considering the systematic uncertainty, an upper 
limit on the product branching fraction is determined to be 
B(v (3686) > yn-(2S)) x B(n(28) + yJ/v») < 9.7x 10-6 
at the 90% C.L., where the systematic uncertainty is incor- 
porated by a factor 1/(1 — osyst.) for conservative. Com- 
bining the result of B(w(3686) — yn-(25S')) obtained by 
BESIII Bh, the upper limit of the branching fraction for 
Ne(28) > «yJ/v is B(n-(28) — yJ /Y) < 0.044 at the 90% 
C.L. Using the width of 7.(2S) of 11.3*2:2 MeV/c? [17], our 
upper limit implies a partial width of I'(7.(2S) > yJ/w) < 
0.50 MeV/c?. Although this result agrees with the prediction 
of LQCD (0.0013 MeV/c?) (38h, it clearly has a very limited 
sensitivity to rigorously test the theory. 
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